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ABSTRACT 



Aims. I present the Monte Carlo radiative transfer code TRADING (Transfer of RAdiation through Dust In Galaxies). The code 
computes self-consistently the extinction of radiation in a dusty medium (including absorption and scattering) and the dust emission. 
Methods. A binary-tree adaptive grid is used for the description of the dust distribution. Dust radiation is computed at thermal 
equilibrium or under stochastic heating condition, for a distribution of grains of different radii and materials. The code is applied to 
the case of a clumpy galactic disk, including both diffuse dust and a distribution of spherical clouds modelled on the GMCs of local 
galaxies. Diffuse and localised sources of starlight are used, with independent spectra. 

Results. A model is provided for the edge-on galaxy NGC 891. The SED of the galaxy from the UV to the submm/mm range can 
be well reproduced by: a bulge/disk configuration of old stars together with an extended dust disk, as suggested by the analysis 
of optical/near-infrared images; a clumpy dust distribution of the same mass as the diffuse dust disk, together with a UV emitting 
component, half of which in the form of a diffuse disk and half in sources embedded in clouds. In total, it is found that about 35% 
of the bolometric radiation is absorbed (and emitted) by dust; and that absorption of starlight from the old population contributes to 
about 60% of the dust emission. A significant component of the dust emission from clouds is due to absorption of diffuse radiation. 
Radial profiles of dust emission in a clumpy disk are almost independent of the wavelength, with the exception of the wavelength 
range on the Wien side of the thermal equilibrium peak. 

Key words, dust, extinction - radiative transfer - methods: numerical - infrared: galaxies - galaxies: structure - galaxies: spiral 



1. Introduction 

The last two decades have seen a good effort in studying the role 
of dust in the appearance of a spiral galaxjQ Radiative trans- 
, fer models have shown that it is necessary to include scatter- 
■ ing and consider appropria te geometries for a proper analysis 
' of t he internal extinction ( Pierini et al. l l2004l:lTuffsetal.ll2004l: 
. iRoc ha et al. 2008, to cite only the most recent works). Building 
on models of extinction in the optical, simulations of dust emis- 
, sion in the Far-Infrared (FIR) have been prod uced (Bianchi et al.l 
. I2000al: iPopescu et alfcOOOt iBaes et al.ll2004l) . 

Puzzling has been the comparison between models and ob- 
servations. Dust disks of moderat e optical depth can r e produce 
obser vations of edge-on spirals jXilouris et al.l 119991; iBianchil 
l2007h . in accordance with the analysis of exti nction of back- 
groun d sources by less inclined disks (see, e.g., iHolwerda et al.l 
12007 ^. Instead, the Spectral Energy Distribution (SED) in 
FIR/submm requires dust absorption and emission from a sig- 
nificantly larger mass (iBianchi et al fcoOOatlPopescu et alj200a 
iMisiriotis et al.ll2001h . 

The excess mass derived fro m emission models co uld be as- 
sociated with molecular clouds dPopescu et al.l 1200 0). Thus, it 
becomes necessary to study the influence of inhomogeneities, or 
clumping, on the radiative transfer. So far, the effects of clump- 
ing have been studied mainly in extinction for galacti c disks 
(Bianchi et al. 2000b; Matthews & Wood 2001; Pieri ni et all 



12001 and in emission for g eometries appropriate to starburst 
galaxies (Mis selt et al.l l200ll) . In all these cases, clumps have 
been modelled as higher density cells in a uniformly spaced dust 
grid. However, the lack of resolution inside a clumpy cell does 
not allow to study cases in which significant emission comes 
from within the cell, as when modelling young stellar sources 
embedded in a molecular cloud. To obviate this, iSilva et all 
( il998.) use a complete radiative transfer solution for a single 
cloud, and add the cumulative output from clouds to the dust 
emission from the diffuse medium. A more self-consistent so- 
lution would require the usage of an adaptive mesh refinement 
grid, with a finer resolution for the denser medium. 

I present in this paper the Monte Carlo (MC) radiative trans- 
fer code TRADING (Transfer of RAdiation through Dust IN 
Galaxie^. For selected viewing directions, the code produces 
images of dust-extinguished starlight in the Ultraviolet (UV), 
Optical and Near Infrared (NIR), and of dust emission from the 
Mid-Infrared (MIR) to the submm. Images can be integrated to 
obtain the total SED. TRADING builds on our previous models 
and has been developed from the regular grid version used in 
Bianchi (2007) . 

The scheme of TRADING is similar to tha t of the code 
DIRTY (iGordon et all 120011: iMisselt et al.1 12001 1). However, in 
TRADING the dust distribution is described with an adap- 
tive grid, while DIRTY is based on a regular grid. An adap- 



' A parallel and more considerable development has occurred in ra- 
diative transfer studies of dusty circumstellar enviro nments (see, e.g., 
the various models described in iPascucci et d]|2004l) . 



- The acronym was suggested by a misprint in a notice at the 
Exploratory Workshop :Tracing Dust In Spiral Galaxies. The workshop 
was funded by the European Science Foundation and took place in 
Ghent, Belgium on May 13-17, 2007. 
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live grid ha s been included in the dust radiative transfer codes 
SUNRISE (iJonssonI |2006|) . which has been used to study 
extinction in hydr odynamical simulati ons o f isolated spirals 
jRocha et al.ll2008! ). and ART- ( Li et al.ll2008l) . SUNRISE does 
not simulate dust emission, while ART^ includes thermal equi- 
librium heating from a single grain. In TRADING, instead, 
emission is fully taken into account by using a distribution of 
grain sizes and materials; it includes both thermal and stochas- 
tic heating, and self-absorption. With DIRTY, SUNRISE and 
ART^ (and many other codes in the radiative transfer literature) 
TRADING shares the use of the MC technique, as this is of easy 
implementation in any geometrical configuration and not nec- 
essarily heavier in comp uting time than analytical techniques 
dBaes & Deionghell200ll) . 

TRADING has been developed, and used in this paper, to 
study the influence of inhomogeneities (clouds) in the dust emis- 
sion from a galactic disk. To my knowledge, this is the first at- 
tempt ever to model the influence of dust on a galactic SED 
by using a fully self consistent radiative transfer calculation 
over scales ranging from molecular clouds to the galactic disk. 
Despite the code is admittedly tuned to this problem, however, it 
is of general use and can be easily adapted to other astrophysical 
scenarios. 

The scheme of the paper is as follows. The code is de- 
scribed in Sect. |2] In Sect. [3] I present the implementation of 
the Draine & Li (2007) dust model used in this paper In Sect.|4] 
TRADING is used to model the global SED of NGC89L The 
results of this application, together with the main features of 
TRADING, are summarised in Sect. |5] 



2. The code 

TRADING consists of two main parts: the first follows the MC 
random walk of stellar energy packets (photons) in the dusty 
medium, through absorptions and scattering; a 3-D map of the 
Interstellar Radiation Field (ISRF) is produced and fed to the 
second part of the code, which computes dust emission. 

The application of the MC technique to radiative transfer 
problems is well described in the literature. Thus, I do not pro- 
vide here full details on all the algorithms used in TRADING 
but I will refer the reader to other published works (in particular, 
to Bianchi et al 1996; Gordon et al. 2001; Misselt e t al. 200lt 
iBaes et al ]|2003tlNiccolini et al.ll2003t lj onsso 3l2006t) . 

2.1. Stellar Emission 

In TRADING, stellar radiation can be emitted from three differ- 
ent kind of distributions: an exponential disk, a spheroidal bulge, 
a collection of discrete spherical or point-like sources. 

As usual in models of spiral galaxies, the luminosity density 
of the stellar disk is given by an exponential distribution along 
the radial, r - -\/x^ + y^, and vertical, z, coordinates. It is 



p^'^Hr,z)=p*^kexp 



kl 



(1) 



where /jj and Zs are the radial and vertical scalelength of the stel- 
lar distribution, and is the central l uminosity density. 

The de-projection of a Sersid d 19681) profile is used for the 
spheroidal bulge. For a profile index n, the luminosity density 
can be written as 



p^-'^^r,z)^p 



buigeexp[/7„(l-B'^")] 



where p'^"'^'' is the luminosity density at B = 1, 
o — , 

/7„ = 2« - 1/3 - 0.009876/«, a^(2n- l)/2«, 

b/a is the minor/m ajor axis ratio and /?e the effective radius 
(iPrugniel & Simien 1997). Bulges with n = 1, 2, 3 or 4 are im- 
plemented in TRADING. For « - 4 the projec tion of Eq. |2]on 
the sky plane is the R^^'^ dde VaucouleursllT959[) surface bright- 
ness profile. 

To simulate emission from star-forming regions within 
GMCs (Sect.|Z2j, TRADING allows radiation to be emitted also 
from a list of spherical sources. Each source is characterised 
by its fractional luminosity relative to the total luminosity of 
all sources in the list, /",; by its radius, r,-; and by its center in 
space. In the present version of TRADING, the luminosity den- 
sity within each source is constant (though I mostly use in this 
paper point-like sources, r,- = 0). This can be easily changed 
to account for a specific radial distribution. For example, the 
SPH smoothing kernel could be used to simulate emission from 
sources in hydrodynamical simulations (Jons son 2006). 

Several emission components can be included in a simula- 
tion. Each component is described by its set of parameters and 
by its fractional luminosity /, with respect to the total model lu- 
minosity. The position of emission of a photon is then derived 
using the MC method. First, a component / is selected so that 



(3) 



k=0 



k=0 



Throughout this paper, is a random number uniformly dis- 
tributed in [0, 1]. Then, the position within the selected compo- 
nent is derived: the procedures to obtain the x,y,z coordinates of 
photon e mission from th e expo nenti al disk and t he bul ge are de- 
scribed in iBianchi et all (Il996h and iBaes et all (l2003h . To save 
computing time, a table with L'"''s'^(r), the luminosity of a Sersic 
bulge emitted within radius r, is provided to the code, and the 
distance from the bulge center is obtained through interpolation. 

If the photon is emitted within the list of spherical sources, 
the specific source first need to be selected. This is achieved by 
applying Eq. [3]to the fractional luminosity of a source The 
x,y,z coordinates are then extracted in a manner analogous to the 
bulge, with the radial distance from the source center given by 
KiR^I^, as results from applying the MC method to a homoge- 
neous sphere. 

The direction of emission of photons is isotropic in the solid 
angle, and the initial photon energy (or weight) is set to unity. 
The photon wavelength A is obtained from the adopted stellar 
spectrum F,i, 



F^dA, 



(4) 



(2) 



with [/^i^i", /Imax] the wavelength range for stellar emission. To 
avoid numerical inversion of Eq. |4] within the code, a table is 
provided for each spectrum. Each of the emission components 
can have an independent stellar spectrum. 

2.2. The dust distribution 

In radiative transfer models of spiral galaxies, dust is usually 
assumed to be distributed in a smooth exponential disk similar 
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to that adopted for stars. At the reference wavelength of the V- 
band, the extinction coefficient of the smooth disk can be written 

as 

J.o. 

r \7\ 

(5) 



kY(r,z) = - — exp 
2 Zd 



r 

hi 



M 

Zd 



where r^" is the central, face-on, optical depth of the dust disk, 
and /id and Zd are the radial and vertical scalelengths. 

If the dust properties are constant within the disk, Eq.|5]can 
be integrated over the volume to give a total dust mass 



Mdust 



G/D TV 



G/D TV 



(6) 



where mu the mass of a hydrogen atom, G/D is the gas (H) to 
dust mass ratio and tv/A^h is the V-band optical depth per unit 
H column density. The las t two quantities are pr ovided by the 
adopted dust model. As in iBianchi et alj ( l2000bl) . I use Ty°' as 
an indicator of the amount of dust in a model, regardless if part 
of the dust is distributed in a clouds. That is, a dust distribution 
labelled by a given value of Ty°- has the same Mdust as a smooth 
model with that central face-on optical depth. 

A fraction of the dust (gas) mass is assumed to be dis- 
tributed in clumps, while the rest is in the smooth disk. Thus, 
in the presence of a clumpy distribution, the smooth disk has 
a central face-on optical depth (1 - /c)t^°-. In TRADING, the 
clumpy distribution consists of a collection of sp herical, homo- 
geneous clouds. Following Bian chi et al.l (l2000bh . their proper- 
ties are modelled on those of Giant Molecular Clouds (GMCs). 

The mass distribution of GMC in the Milky Way and in a few 
other Local Group galaxies is well described by a single power 
law, 



(7) 



oc M ^ , 

dM 

over a wide range of cloud masses (iBlitz & Wilhamsl [1991 
iBlitzet al.ll2007l) . For y < 2, the H mass of a cloud / can be 
derived with the MC method from 

with Minf and M^^^ the minimum and maximum mass allowed 
in the spectrum. The total number of clouds A^c is determined 
by requiring the conservation of the dust mass in the clumpy 
component, i.e. 

(G/D)-'YjMi« fcMdu^ 



(due to the random nature of the cloud mass derivation, exact 
equality is not assured). 

Observations also suggest that the gas surface density of the 
clouds is independent of the cloud mass, with a mean value E 
that changes by less than a factor of two f rom galaxy to galaxy 
jBlitz et al.ll2007t iRosolowskv et al.ll2003l) . If a constant 2 is as- 
sumed, the radius R follows directly from the cloud mass. A con- 
stant 2 also implies a constant optical depth along a cloud diam- 
eter. For a homogeneous spherical cloud, the V-band extinction 
coefficient is then 



3 S TV 



(8) 



4 ninR Nh 

As for the clouds spatial distribution, it is assumed that it 
follows a double exponential as in Eq. |5] The location of each 
cloud is also used to define the location of the spherical sources 
of radiation described in Sect. 12.11 



2.3. The dust adaptive grid and the random walk 

The analytical description of the dust distribution outhned so far 
is used to fill up the dust adaptive grid. A 3-D binary tree (an 
octree) is constructed subdividing recursively each volume el- 
ement (a cell) into eight children. First, the dust distribution is 
enclosed in a cubic volume of size D (the root cell). Then, the 
condition 



f kydV<E I 

Jcell Jd' 



k^/dV, 



(9) 



is checked, with £ < 1 a grid re finement parameter dWolf et al.l 
ll999l:lKurosawa & Hillieill200lh . I have used E = 10"^^ for the 
simulations presented in this paper. If the condition is not met, 
the cell is split into eight cubes and the condition checked again 
recursively until it is met. If the condition is met, subdivision 
is ended and a leaf cell is found. If L subdivisions have taken 
place, the tea/ cell has size D/(2^)0. For each cell, information 
on its boundaries are stored, together with pointers to its father 
cell (unless it is the root cell) and to its eight children (unless it 
is a leaf cell). The local value of fcv is assigned to each leaf cell. 

The criterion in Eq. |9]is based on the mean opacity of a cell 
(equivalently, on the density, if the dust properties are identical 
throughout space, as in this work). Since it is not tied to the local 
radiation field (which is not known before the radiative transfer 
calculations), it may fail in producing an adequate cell resolution 
in the presence of very strong field gradients. A high resolution 
is needed for the case of clouds with internal point sources, and 
could in principle be achieved by using a smaller value for E. 
However, this would result in a larger number of cells and in 
an unnecessary finer sampling of the diffuse mediunfl To avoid 
this, the grid algorithm allows to force splitting in cells: if the 
cloud diameter is smaller than the dimension of a cell of level L, 
the leaf cells can be required to split further by I sublevels, down 
to level L + I, thus providing at least 2'"' resolution elements 
along a cloud diameter. An example of the binary structure of 
the dust grid is shown in Fig.[T] 

The optical depth of the dust distribution along any given 
direction is found from the intercepts of the photon path with 
the boundaries of the leaf ce lls in the grid (for more details, 
see iKurosawa & Hillierl I200TI) . The leaf cell that contains the 
position of emission is found by recursively ascending the 
bin ary three from the root: the point location algorithm of 
iFrisken & Perrvl (|2002|) is used. The coordinates of the intercept 
of the photon path with the boundaries of this first cell are used 
to locate the next cell along the path. This neighbour cell could 
be found again by ascending the tree from the root. An alterna- 
tive scheme would consist in finding the smallest non-Zea/cell 
that is a common ancestor for both a cell and its neighbour: the 
tree is then descended from the first cell to the common ancestor, 
and then ascended again from the ancestor to the cell neighbor. 
This second implementation is more efficient (Frisken & PeiT^ 



^ As an example, let's consider a homogeneous root cube: setting, 
e.g., E = 8"^ is equivalent to produce a regular grid of 8' fea/ cells (de- 
scending down to level L = 3), each of size D/8. Using a cubic root vol- 
ume to enclose a thin disk-like structure may seem inadequate, since re- 
gions distant from the galactic plane are empty. However, these regions 
are enclosed by a limited number of low-level leaf cells. Furthermore, 
the cubic geometry makes the structure of more general application and 
is more effective in the description of the spherical clouds in the clumpy 
distribution than cells with a reduced size along z- 

For smooth models, and for the smooth medium in the clumpy mod- 
els presented here, the chosen E value ensures enough spatial resolution 
for the convergence of the results. 
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Fig. 1. A cut through the binary dust grid used in Sect. 14.21 per- 
pendicular to the galactic disk. Since mean values are stored 
in each cell, the density on the border of clouds does not drop 
abruptly to the value of the smooth medium. The size of the box 
is 1.5 kpc X 1.5 kpc. 



l2Q02 h. If dl is the length of the intercept, the optical depth along 
the path in the cell is 

where A is the wavelength of the photon and is the dust extinc- 
tion law from the adopted model. The process is then repeated 
and the optical depth along the path accumulated, until the 
photon suffers the first scattering event (or the boundary of the 
root cell is reached). 

The location of the scattering event is determined randomly, 
by finding the optical depth 

T=-log(l-!R), (10) 

and stopping the path in the grid where = t. If t is larger than 
the total optical depth along the chosen direction, the photon es- 
capes dust. To ensure that all photons contribute to the scattered 
flux, e ven in the c ase of low global opacity, scattering may be 
forced dCashwell & Everett 1959) . Typically, the first scattering 
event is fo rced, while the foUowing are determine d according 
to Eq. Ql ("B ianchi et al] 1 19961: iGordon et alj|200ll) . The same 
scheme is adopted here. 

In the scattering position, the photon propagates along a new 
direction, which is derived rando mly from the adopted s catter- 
ing phase function. Here I use the lHenvev & GreensteinI (Il941h 
phase function, dependent on the scattering asymmetry parame- 
ter gx - (cos(0)), with 6 t he angular deviation be tween the new 
and the old directions. The lHenvev & GreensteinI phase function 
is a reaso nably good app roximation for the wavelengths consid- 
ered here (lDrainel2003bl) . The values for coa and gA are provided 
by the adopted dust model. 

The scattered photon has its weight reduced by the albedo 
loa, while the fraction (1 - (Oa) is absorbed. The absorbed en- 
ergy could be stored in the /ea/cell where scattering takes place 



and used later to derive dust emission dBianchi et alj l2000al: 
iGordon et al.ll200TI ). However, it is more efficient to store the 
absorbed energy along all the cells crossed by a photon rather 
than at the end po int only, a con c ept in troduced bv lLucvl d 19991) 
and described in Niccolini et al.l (l2003i) for both the non-forced 
and forced scattering cases. Using this formalism, TRADING 
accumulates the amount of energy absorbed as a function of the 
wavelength, Wa, from which the local ISRF is derived, 

1 

AnkAd-cDA) AV' 

with Ay the cell volume. For each leaf cell, the average of J a 
is stored in A^ggp contiguous wavelength bins, logarithmically 
spaced in the range [A^^\ A^^H 

After the first scattering, the whole procedure is reiterated, 
until the photon weight falls below a limit value (lO"'* is used 
here). A number A^phot of photons are run to provide high signal- 
to-noise images and Ja values in the grid. The model can be 
normalised to physical values by multiplying the outputs by 
Lbo\/Npho(, Lbo\ being the bolometric luminosity assumed for the 
stellar sources. 

Simulated images of dust-extinguished starlight can be also 
produced, for all the N^^^ bins or at specific values of A. Any 
viewing directio n is possible. This i s achieved using the peeling- 

technique of lYusef-Zadeh et aL l (IT984l) . in which each emis- 
sion and scattering event contributes to the surface brightness 
distribution. The technique allows to produce high signal-to- 
noise results without the need of collectin g photons in a broad 
angle band (see, e.g.. iBianchi et al.lll996l) . thus preventing im- 
age smearing when there is a strong gradient with the viewing 
direction (as for, e.g., a disk seen edge-on). 

Once that the radiation field has been determined, it may re- 
sult that some cells contribute negligibly to the total absorbed 
(and emitted) energy. To avoid unecessary calculation, the grid 
is reduced in size by adopting a criterion similar to that in 
iMisselt et al.l (12001 1) : cells are ordered according to the amount 
of energy absorbed within them, and a threshold is derived so 
that all cells below the threshold contributes to a fraction /abs 
of the total absorbed energy. The grid is recursively accessed 
to find fathers of leaf cells which have all eight children below 
the threshold. When such cells are found, the children are cut 
and their amount of absorbed energy is transferred to the father, 
which becomes leaf. Instead, when leaf cells below the thresh- 
old are isolated, they are excluded from the calculation of dust 
emission. When a clumpy dust distribution is considered, even 
a small value of f\,s may result in an appreciable reduction of 
the memory occupation of the grid and of execution time of the 
emission code. In this work, I have used /abs - 0.1%. No sig- 
nificant change in the results is produced by the adoption of this 
threshold. 

2.4. Dust emission 

When considering the transfer of radiation through absorption 
and scattering, it is not necessary to know the details on the dust 
composition and size distribution, but it is sufficient to consider 
mean properties (Wolf 2003) : these are the values Aa, a>A and gA 
used in the previous section. 

Dust emission, instead, depends on the properties of each 
grain. For an ensemble of spherical grains of A^^at different ma- 
terials, each characterised by a distribution nt(fl) for the radius a, 
the emission coefficient (energy emitted per unit volume, time. 
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wavelength and solid angle) for thermal equilibrium emission 
can be written as 



k=l ^ 



Ukia) a' QT{a,A) BM da 



"mm p 

^ I ntia)a^ Qf{a,y) da. 



(11) 



k=\ 



with 2^"' and Qf'^ the extinction and absorption (emission) ef- 
ficiencies (i.e. the cross-sections divided by na^), Bi the Planck 
function and the equilibrium temperature of a grain, depen- 
dent on the radius a and on the material k. In TRADING, the 
size distribution for each material k is sampled with A^,^ contigu- 
ous bins logarithmically spaced between the minimum and max- 
imum radii of the adopted grain distribution. Equation [TT] thus 
becomes 



A'mal K 



k=l j=i 



3jk^BATi(aj,k)), (12) 



where Qj^ is the absorption cross section integrated over the j- 
th radius bin of width Aaj, normalised by the V-band extinction 
cross sectior|3. 



Qjk- 



rikiaj) a J Ql \aj,A) Aaj 



(13) 



k=i ^ 



Ukia) a^ QT^ia,^) da. 



At thermal equilibrium, can be found equating the radia- 
tion absorbed from the ISRF J,\ to the dust emission. For the y-th 
grain of material k, the thermal balance can be written as 

Y^J;^Qjk{A)^Ai^ j BMaj,k))Qjt(A)dA, (14) 

where A/l, is the width of the i-th bin used to store J^. To save 
computing time, the integral in Eq.[T4]is tabulated, and the tem- 
perature Tiiflj, k) is found by interpolation. 

For small grains (a ^ O.Olyum) heating is stochastic and it 
is not possible to define an equilibrium temperature. Instead, a 
given grain is characterised by a temperature distribution PiT), 
dependent on the grain c haracteristics and on the heating field. 
Following the method of iGuhathakurta & Drain"^ ( 17989,) . P(T) 
can be derived from the transition matrix A fj, that gives the prob- 
ability per unit time that a dust grain is heated (or cooled) from 
the enthalpy state //, to the enthalpy state ///.In the notation 
used so far, for a given grain, the matrix elements due to heating 
by the ISRF are (/ > /) 



Qjk(A) AHf 
(Hf-Hiy ' 



where AHj is the width of the / enthalpy state, A - hc/(Hf-Hi) 
is the wavelength associated to the transition, and the value of 
is obtained by interpolating the ISRF provided by the radiative 
transfer code (see Misselt et al. 2001, for a full description of 
the matrix elements). Considering only cooling terms from level 



^ A table with QjkW is computed from the adopted dust model and 
provided to the radiative transfer code. Ancillary tables are produced 
by averaging Q over the wavelength bins used for the sampling of j.i 
(Eq.IEll, J A (Eq.[Tl and 7f -. 



f+ 1 to level / (th e thermal continuous cooling approximation; 
IDraine & LillMl . it is 



/ QjkW BATf^i)dA, 



f 



where Tf+i is the temperature associated to state f+l. In 
TRADING, a temperatures range for P(7'd) is selected, and A^j- 
states are defined, logarithmically spaced in temperature. The 
enthalpy associated to each state is derived by integrating over 
the temperature the specific heat adopted for the grain materials. 
The same temperature range and Nt are used for all grain, re- 
gardless of the local ISRF intensity and spectrum. This allows to 
compute tables for the cooling terms and for the multiplicator of 
in the heating terms beforehand, thus saving computing time. 
For a grain of material k and radius a y that undergoes stochastic 
heating, the Sjk term in Eq. [T2]becomes 



Nt 



^Jk 



(15) 



In each cell, the emission coefficient is computed by sum- 
ming the contribution of all sizes and materials in the dust distri- 
bution. The calculation is carried out for contiguous wave- 
length bins, logarithmically spaced in the range [A'^^^, A'^^}. 

Subsequent calculations are carried out using a no-scattering 
MC radiative transfer procedure. The MC allows to take into 
account self-absorption, i.e. the possibility that dust-emitted ra- 
diation is again absorbed by dust. This could be relevant in some 
astrophysical situations, when MIR radiation travels through a 
very dense dusty medium. 

In principle, the cell where a photon (of unit weight) is emit- 
ted should be selected randomly from the amount of energy 
emitted within it, and the wavelength determined randomly from 
the local spectrum. However, it is not straightforward (or effi- 
cient in terms of computing time) to access the binary tree in this 
way. Furthermore, the procedure would also require to store i,\ in 
the memory for each cell. Also, since scattering is not included, 
calculations can be carried out at once for all wavelengths. Thus 
in TRADING, once is derived, a few polychromatic photons 
are emitted (five are sufficient for the simulations presented in 
this paper), with a /(-dependent weight proportional to jaAV {AV 
is the cell volume). The emission coordinates within a cubic 
(homogeneous) leaf cell and the direction propagation are de- 
termined in the usual way. As the photon traverses cells, a local 
ISRF JY' can be computed from self-absorption in the wave- 
length range /Imaxl' analogously to what done for stellar 
photons. In case, the wavelength range could be limited by a 
maximum wavelength A\:^^^, if self-absorption is insignificant for 
A > /i^ax' and a different number of bins over which is 
stored, A^ggp, can be selected. Selecting a restricted wavelength 
range for self-absorption may become necessary when computer 
memory is a concern, since, contrary to j^, it is necessary to store 
ys.a. fQj. ^jj ^gjjg before carrying on the calculations. 

Once " is derived, its contribution to dust heating is added 
to that of J^, and a new emission coefficient j,i can be computed. 
The entire process is repeated until the total output of dust emis- 
sion (i.e. the dusty SED integrated over the solid angle) does not 
change by more than a given fraction, typically a few percents. 
The convergence check, as well as the splitting of TRADING in 
the extinction and emission parts, could be avoided by using the 
immediate re-emission concept of Biorkman & Wood (2001'). 
With this procedure, the absorption of a stellar photon is im- 
mediately followed by the emission of an infrared photon from 
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Fig. 2. Tests of the T RADING code. On the left, comparison with the DUSTY solutions of the benchmark problems of 
llvezic & Elitzui (|1997 ) for a hom ogeneous sphere. On the right, comparison with the RADICAL solutions for the benchmark 
problems of Pascucci et al.l (|2004|) . In all cases, the SED has been normalised to the total luminosity; the relative difference has been 
computed with respect to the benchmark solutions. See text for details. 



dust, with the wavelength sampled from a distribution that cor- 
rects itself at each absorption event; eventually the distribution 
of infrared photons converges to the correct distribution, without 
the need for a separate calculation of the absorption and emis- 
sion processes. Although the method has a solid mathematical 
ground in the case of thermal equilibrium heating and in MC 
procedures where phot ons have wei ghts, it is not clear yet if it 
can be applied with t he iLucvl jl^ Qg*) estimato r for 7^ o r in case 
of stochastic heating (iBaes et al.l ll2005; but see lJonsso for 
a different view on the second issue). Because of this, I have 
preferred the current scheme. 

As for starlight, images of dust-emitted radiation can be ob- 
tained from any viewing direction, for all the Nf^^ bins or at 
specific values of A, using the peeling-off technique from each 
emission point in the MC procedure. 

2.5. Benchmarking 

The code results have been tested on a few radiative trans- 
fer solutions^_die outputs from the MC radiative transfer code 
of iBianchi et al.l (fl996) : the description of the radiation field 
due to a point source in th e center of a finite scattering sphere 
dSiewert & Maiorino' 1979); the benchmark problems in ID ge- 
ometry of Ivezic et aL (1992); th e benchmark problems in 2D 
geometry of iPascucci et al.l (|2004|) . I describe here the compari- 
son wdth_flieJa^ttwobenchmark cases. 

llvezic etal] ( Il997h describe a set of radiative transfer prob- 
lems for a spherical dust distribution heated by a central point 
source. I have compared the output of TRADING with the case 
of a homogeneous sphere. The solutions have been obtained in 
numerical format using the publicly av ailable radia tive transfer 
code DUST'ifl based on the method of llvezic & EUt zur ( 1997^ 
DUST Y is one of the radiative transfer codes used in llvezic et al.l 
(IT997h . 

^ |http : //www ■ pa . uky . edu/"moshe/dusty/| 



A grid for a homogeneous sphere of radius Rout was first con- 
structed. The dust distribution of the benchmark cases has a inner 
cavity for Rin - /?out/1000. Most of the sphere volume is sam- 
pled with cells of level L - 7 (thus, having a size of 2 xRo^t / 1 28), 
but finer subdivisions (down to L = 14) are used for a proper 
rendering of the inner and outer boundaries. In total, the grid 
consists of 6 X 10^ leaf cells. The extinction coefficient (at the 
reference wavelength of 1 yum) is given by the optical depth of 
the model along the radius, ki = Ti/iRo^t - Rin)- 

DUSTY only deals with a single grain of mean properties 
and the scattering is assumed to be isotropic. Thus, = has 
been set in TRADING , and and co^ have been chosen as in 
llvezic &Ehtzu3(II997l) . 

In the benchmark cases, radiation comes from a central point 
source, emitting as a blackbody at - 2500 K. That ISRF has 
been sampled with 50 bi ns in the range [0.3,2 0] fim. DUSTY 
solutions are scale-free (llvezic &Elitzuj[T997i) . and it is not 
necessary to provide the absolute luminosity. Instead, for the 
comparison with TRADING, a physical value for i?out was as- 
sumed, and the bolom etric luminosity Lhoi wa s deri ved with the 
relatio ns provided in llvezic & Elitzuj ( Il997l) and llvezic et all 
( Il997l) . Dust radiation has been sampled with 100 bins in the 
range [1, 1000] fim. The same sampling has been used for the 
contribution of self-absorption to the ISRF. A total of 10^ pho- 
tons have been used. 

The result from the comparison is shown in the left panel 
of Fig. |2] for the cases t, = 1, 10 and 100. As TRADING 
deals with extinction and emission separately, the stellar and dust 
SEDs have been summed together in the overlapping wavelength 
range. The inclusion of self-absorption has been necessary for 
a correct solution of the two high optical depth cases: conver- 
gence to 0.5% was required. The agreement with the benchmark 
cases is satisfactory, with TRADING differing from DUSTY by 
only a few percents in most cases. The agreement is worse at 
shorter wavelengths, for several reasons: the intrinsic noise in 
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the MC procedure; the difficulty of properly sampling the inte- 
gral of a steeply rising spectrum (Eq.|4| with tables; the lack of 
ad-ho c optimisations for high optical depths (see, e.g., Juve^ 
l2005h . In particular, TRADING diverges significantly from the 
DUSTY solution for the Ti = lOOcase, butonly for/l S 3/vmand 
very low flux levels 10"^ in the units of Fig.|2]i. Calculations 
for these low fluxes in the higher optical depth case also pose 
problems for DUSTY: numerical errors had to be prevented by 
requiring a significantly higher numerical accuracy than in the 
other cases. Nevertheless, the radiation field appears to be com- 
puted correcly by TRADING in all cases, as shown by the bulk 
of dust emission in the SED. Also, the surface brightness profiles 
of stellar and dust radiation show a similar agreement. 

The geometry analysed so far is equivalent to what will be 
used in Sect. 14.21 for embedded stellar sources at the center of 
GMCs. Unfortunately, in simulations of galactic disks it is not 
possible to have, for each cloud, the same level of spatial reso- 
lution as in the comparison presented here. Thus, it will become 
necessary to check for the convergence of the simulation outputs. 

A benchmark problem with a more compelling geometry is 
provided by Pascucci et al. (2004). They use an axially symmet- 
ric dust disk, extending from an inner cavity of 1 AU to an exter- 
nal boundary of 1000 AU, with a steeply rising density gradient 
close to the center The dust grid I have used for this problem has 
been obtained with the procedure described in Sect. 12. 31 but forc- 
ing cell splitting in the inner 50 AU down to L = 14 close to the 
center. The grid consists of 1.1 x 10^ leaf cells. Dust properties 
are provided by the authors, and isotropic scattering is assumed. 

As for the previous case, the heating source is point like and 
located in the center, with Lboi = l^^o and blackbody emission 
at Ti, = 5800 K. For this case, the ISRF has been sampled with 
50 bins in the range [0.12, 10] fim. All the other parameters are 
the same as for the comparison with DUSTY. 

As a reference, I have used the solution of the benchmar k 
problems provided by R ADICAL (|Du llemond~ & Turollall2000l) . 
one of the codes used in lPascucci et al.i (i2004i) . The result of the 
comparison with their edge-on case (inclination / - 77.5° from 
the polar axis) is shown in the right panel of Fig. |2l Again, the 
agreement is satisfactory. For most wavelengths the TRADING 
solution is within 10% of the RADICAL solution, a difference 
only sligthly larger than wh at achieved by compar ing the outputs 
of the various codes used in lPascucci et alj (|2004|) . Only the high 
optical depth model shows a larger discrepancy on the stellar 
side of the SED. This is due to the limited resolution of the grid. 
Better agreement is o btained for the more face-on models. 

As pointed out bv'J onssonI (l2006l) . regular or adaptive 3D or- 
thogonal grids are not ideally suited for a comparison with the 
benchmark solutions described here. In fact, no benchmark so- 
lution is available yet for the case of diffuse sources cospatial 
with dust, a configuration more suitable for the study of internal 
extinction and emission in galaxies. This is the subject of a re- 
cently started project on comparing different codes that cater to 
the study of radiative transfer in galaxies. 



3. The dust model 

The dust grain model of iDraine & Lil (|2007|) has been adopted 
for the simulations presented in this work. The model consists in 
a mixture of silicate and carbonaceous grains. The size distribu- 
tions was derived by fitting the extinction law and infrared emis- 
sion from dust as observed in our Galaxy (Weingartner & Draine 
l2001al;lLi & Drainel200ll;lDrainel2003al flie model for the Ry = 
3.1 extinction law is used here). 



The dust grain extinction and absorption cross-sections have 
been derived using the Mie theory for spherical particles (Mi^ 
119081 ; see also Bohren & Huffman 1983). The optical proper- 
ties of amorphous astronomical silicates are used for the sili- 
cate component and those of graphite for the larger carbonaceous 
grains. 

For small carbonaceous grains, the absorption cross sec- 
tions are derived from the properties of PAH materi als, the most 
probable responsible of the MIR emission features dLi & Draind 
rame & Lill2007 ). In the model used here, these grains 
contain 4.6% of the mass of the carbonaceous components. The 
properties of PAH grains als o depend on the ionization state, 
and the ionization fraction of iDraine & Lil (12007) is used. The 
absorption cross section of the carbonaceous components is as- 
sumed to change smoothly from that of PAHs to that of graphite 
at a radius of 0.005 fim: the smooth transition provides a contin- 
uum absorption component to the PAH features for small grains, 
whil e for a ^ 0.01 nm the graphite component becomes domi- 
nant jDraine & Lil l2007l) . 

With the above recipe, I have computed the mean extinction 
law A A, albedo w,i and asymmetry parameter g,i which are used 
for the radiative transfer of stellar radiation (A^ is also used for 
dust radiation if self-absorption is included); and the elements of 
the Qjit table of Eq.[T3] which are necessary to compute the dust 
emission spectrum. In the table, A^mat = 4 grain components are 
considered: silicate grains, with A^^'' = 18 bins in radius from 
a = 3.5A to 0.5;um; graphite grains, with A^,f ^ = 10 bins in 
radius from a = O.Olyum to lyum; ionised carbonaceous grains 
(smoothly changing from PAHs to graphite), with A^^'^^ - 9 bins 
in radius from a = 3.5A to O.Olyum; neutral carbonaceous grains 
with the same bins as for ionised grains. The values of A^* results 
in similar radial bins for both carbonaceous and silicate grains. 
The sampling of the different components was chosen by means 
of trial and was dictated by the necessity of having the smallest 
possible number of bins (since the computing time depends on 
the sums in Eq. [TST i still able to provide a reasonably accurate 
estimate of dust emission. 

Computing time is an issue especi ally for grain s that u nder- 
goes stochastic heating. As noted by Mis selt et al.l ('2001), this 
is the real bottleneck of the computation. For the models pre- 
sented here, the calculation of the temperature probability dis- 
tribution P is limited to grains with a < O.Ol/zm (26 bins out 
of 46), while the rest are assumed to emit at thermal equilib- 
rium. The enthalpy states are sampled wit h A^j- = 80 temper- 
ature bins, with 2.7K< T <2000K. The IDraine & Lil (|200TI) 
models for the specific heats of bulk silicate and graphite have 
been used. While this temperature sampling is sufficient for the 
smaller grains, it does not provide correct solutions for the larger 
grains exposed to stronger ISRFs, whose P is close to a delta 
function. When this happens, thermal equilibrium emission is 
used. A more accurate calculation would require a different tem- 
perature sampling tailored over the heating conditions in each 
cell and for e a ch gra in, as in the iterative algorithm described by 
iMisselt et al.l (12001 1) . However, this is unpractical for the large 
grid considered here, as the method prevents the usage of pre- 
computed tables for the terms in the transition matrix and would 
result in too long computation times. 

Despite the approximations and the limited grain sampling, 
the method used in TRADING results in a reasonably good de- 
scription of dust emission. Fig.|3]shows the emission when heat- 
ing is du e to the average Galactic I SRF (Math is et al. 1983; see 
also We ingartner & Draind 1200 fbl) . scaled by different factors 
U. In the wavelength range where emission is due to stochasti- 
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Fig. 3. Dust emission spectra for the dust model of Sect. [3] nor- 
malised to the hydrogen column density. The heating spectrum 
is the average Galactic ISRF, with intensity scaled by factor U. 



cally heated grains, the agreement with the analogous models of 
iDraine & Lil (120071) is within ^ 10%. A similar agreement with 
full models is achieved using the spectra of Sect. |4] 

Finally, the dust model is characterised hy G/D = 90 and 
tv/Nh = 4.9 X 10-^^ cm^. 




X [ |im ] 

Fig. 4. The stellar SEP fro m the Sb template of 
Fioc & Rocca-Volmerangd ( Il997h . normalised to the bolo- 
metric luminosity. The total SED is shown together with the 
separate contribution of stars of low (M* < 3 Mq) and high 
(M* > 3 M0) mass. For the sake of presentation, the spectra 
have been smoothed over the wavelength. The average Galactic 
ISRF is also shown for comparison. 

1.3 mm (iGuelin et al.lll993h . These fluxes are shown as data- 
points in Fig.|5] 



4. An application: NGC891 



4. 1 . The SED of smooth disks 



In this section, I present TRADING models for a galactic disk. 
The physical quantities assumed in the simulations are tied to 
the case of NGC891, an Sb-type edge-on galaxy whose SED 
and surface brightness distribution s have been subject of ex- 
tensive stud ies (Xilour is^etal.lll999HPop"escu et al.ll2000ll2004l: 
lAlton et ani2b04i: .Dasyra et al.ll2005T) . In accordance with these 
works, a distance of 9.5 Mpc is assumed. 

For stars, I have adopted a typical spectrum of an evolved 
galaxy. It has been ob tained by running the spectral sy nthe- 
sis code PEGASE.l of iFioc & Rocca-Volmerangd (1 19971) with 



the p arameters for their Sb template (see also iLeitherer et al.l 
Il996 ). but without including a dust correction. The intrinsic stel- 
lar spectrum is shown in Fig.|4l it is similar in shape to the aver- 
age Galactic ISRF. By setting an arbitrary cut at 3Mq, the spec- 
trum has been decomposed into the separate contribution of low 
mass stars (emitting most of 0.3//m radiation, and contribut- 
ing to about 80% of the bolometric luminosity), and high mass 
stars (emitting most of UV radiation). This will allow to use sep- 
arate spatial distributions for young OB stars and for the older 
population. 

In all the models, the stellar is sampled with A^^gQ = 20 
bins in the range [0.09, Slyum. Instead, dust emission is sampled 
with A^jgQ = 150 bins in the range [2, 1500]/im. For all the wave- 
length bins, images are produced for selected inclinations, and a 
global SED is obtained by integrating over the maps. 

Models are compared to the SED of NGC891, which is 
known ove r a broad wavelength r ange: tot al fluxes are available 
in the UV jGil de Paz et al.ll2(jo7l). op tical jde Vaucouleurs et al.1 
1199 lb and NIR (iJarrett et al.ll2003l) : du st emission has bee n 
observed at 12, 25, 60 and 100 //m dSanders et al.1 l2003h . 
170 and 200 ^m dPopescu et al.1 12004|) . 260, 360 a nd 580 fim 
jDupac et al.l l2003l) . 450 and 850 yum (lAlton et al.1 [T99&) and 



The structural properties of NGC 891 have been determined 
in the optical/NIR bands by IXilouris et al.l ([1999) (see also 
iXilouris et a 01998 '). Surface brightness fits suggest a thin stel- 
lar disk with h^/zs ~ 10, and a radial scalelength increasing 
from the K to th e B-band (a w ell known feature in less inclined 
object; see, e.g.. lde Jondll996l) . For simplicity, I neglect the gra- 
dient and assume — A kpc, a value appropriate for the NIR 
bands, where the bulk of stellar emission is; Zs = 0.4 kpc is used 
for the vertical scalelength. 

The fits of IXilouris et al.l (1 19991) also provide a description 
of the dust distribution. Since the radial scalelength of the dust 
disk almost doubles the NIR stellar scalelength, I have assumed 
hi - 2/is. The vertical scalelength is taken to be za - Zs/2, and 
Ty° =1 as suggested by the same fits. Similar properties are 
shown by dust dis ks in other edge-on galaxies dXilouris et al.l 
|1999t|Bianchr2007). 

The modelled SED are shown in Fig. |5] (left panel). The 
stellar radiation is assumed to have a bolometric luminosity 
Lboi = 6 X 1O''^L0, so as to provide a match of the optical/NIR 
observations. The result for a single stellar disk is shown for 
a full calculation including stochastic heating (dot-dashed line) 
and for a thermal equilibrium calculation using a single dust 
grain of mean absorption cross section (dashed line). While the 
simpler calculation is not able to produce significant radiation at 
A < 40//m, the output at A ^ lOO/zm is consistent with that of 
the more complex case, clearly showing that most of radiation 
in the dust emission peak is due to thermal equilibrium heating. 
It is necessary, however, to exclude from the equilibrium calcu- 
lations the amount of energy (30% for the case shown in Fig.lSj 
which is absorbed by smaller grains ( Bianchi et al. 2000a). In 
the following, all calculations will be done for the complete dust 
model of Sect. [3] 
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Fig. 5. SEDs for edge-on models. The left panel shows the models for a Ty° = 1 and Lboi = 6 x 1O'"L0. As TRADING computes 
separately dust extinction and emission, the total SEDs have been obtained by interpolating the two outputs over the common 
wavelength range and summing them together. Dust emission only is shown for the case of thermal equilibrium emission by a single 
grain with mean properties (dashed line), after correcting for the energy absorbed by small grains. The right panel shows the SEDs 
of = 2 models. In all cases, the luminosity emitted by low mass stars is 6 x 1O'"L0. In both panels, data points represent the 
observed SED of NGC891 from literature (see text for references). UV, optical and NIR fluxes have been corrected for foreground 
Galactic Extinction. When not shown, error bars are smaller than symbol size. 



An important contribution to ISRF and to dust heating 
may come from a galactic bulge. Following the analysis of 
IXilouris et al.l (Il999l) . I have included a Sersic n - 4 bulge with 
b/a = 0.6 and Re - 1 kpc, with a luminosity 0.3Lboi- The SED 
for this model is shown in Fig. |5] (solid line in left panel). The 
central concentration of the bulge results in an increased tem- 
perature for the larger grains: while in a single stellar disk model 
the larger carbonaceous grains attain a maximum temperature of 
about 28K in the center (25 K for silicates), in the model includ- 
ing a bulge the central regions have graphite grains at 65 K (and 
silicates at 53 K). This is evident in the shift of peak of thermal 
emission towards shorter wavelengths. Changes are instead neg- 
ligible for the stellar SED. The bulge/disk configuration will be 
used in the rest of the paper 



The ■ = 1 disk shown so far grossly underestimate 
(by almost a factor two) the FIR output, as it is able to ab- 
sorb (and reemit) only about 18% of the total energy. The dis- 
crepancy is the well known energy balance problem: a mod- 
erate optically thick dust disk as suggested by optical obser- 
vations is not able to absorb enough energy to match the ob- 
served infrared outpu t; absorption from a more opaque dust 
component is needed ( Bianchi et al. 20004 iPopescu et alj|2000l; 
iMisiriotis etal., 2001) . A = 2 disk with Lboi ~ 7.5 X 10^° L© 
(not shown) can reproduce better the dust emission, but is still 
20% less bright than observations around the peak. However, 
such high optical depths are not not c ompatible with observa- 
tions (IXilouris et alJll999l:lBianchill2007h . 



The SEDs also fail in reproducing the MIR radiation and 
show an excess of UV radiation. Although the second problem 
can be due to an overestimate of the UV component in the in- 
trinsic SED, it is more likely the result of an underestimate of 
opacity: smooth distributions for dust and stars cannot easily re- 
produce the localised extinction of young objects close to their 
dense, more opaque, birth environment. 



4.2. The SED of clumpy disks 

Since the r^" = 1 disk of Xilouris et al.' (1999') fails in repro- 
ducing the FIR SED, Popescu et aL, (2000) have introduced a 
second dust component, of mass similar to that of the first disk, 
but with a radial distribution equal to the stellar This second, 
smooth, disk is devised to mimic the extra extinction in molec- 
ular clouds, while the first could be thought as associated to the 
atomic gas. Indeed, the Interstellar Medium (ISM) in NGC891 
is made by a broad H i distributio n and a more cent rally peaked 
molecular component (see, e.g., lAlton et ani2000h . The distri- 
bution of molecular gas is generally found to be close to the 
stellar (Rega n et al.l 1200 Ih . The molecular and atomic gas in 
NGC891 have similar ma sses (in total 8 x lO'^MfT), 70% of which 
residing in a thin disk; ISofue & Nakail ll993t10osterloo et al.l 
'2007), in agreement also with the mean properties of Sb galaxies 
(Young &Scoville 1991). 

In this paper, the extra mass component required to solve the 
energy balance problem is added in the form of a distribution of 
clouds. The clouds are distributed in an exponential disk of radial 
scalelength h^. The vertical scalelength is the same as for diffuse 
dust. The clumpy component has the same mass as the diffuse 
dust (/c = 0.5). The total dust mass in the grid is Mdust - 1O**M0. 
For the assumed G/D, the dust mass is roughly consistent with 
the gas mass derived from observations. 

In the model, clouds are spherical, homogeneous and their 
opacity is computed from the gas surface density observed in 
GMC. Asssumi ng E - 100 M(T| pc~^, a typical value for L ocal 
Group galaxies (iBlitz et al.ll2007l iRosolowskv et al.ll2003h . the 
optical depth along the diameter is t^'™ ^ 9.2, regardless of the 
cloud dimensions (as follows f rom Eq . [8] using the dust prop- 
erties of Sect. [3]l. According to iRosolowskvl (l2005b . the largest 
clouds have a hydrogen mass Msup = 10^-^ Mq, which corre- 
spond to a radius R ^ lOOpc. GMCs in the Milky Way and a 
few other Local Group galaxies show a power law distribution, 
descri bed by Eq. [7]wit h y = 1 .7 (with the notable exception of 
M33; iBlitz et al.ll2007l) . A value of y < 2 ensures that most of 
the molecular mass resides in high mass clouds. 
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Memory requirements may become heavy for the models of 
a clumpy disk, when a reasonably large resolution is required for 
a large number of cloudfl In particular, a lower mass limit Minf 
implies a larger number of lower mass clouds. Several tests were 
carried out with different values of Mi„f, spanning up to two or- 
ders of magnitude in mass. The results were found to be almost 
independent of the choice on Minf, because of the constant op- 
tical depth and of the luminosity of embedded radiation, which 
was chosen to scale with the mass of the cloud. Thus, a single 
cloud mass is assumed (Minf - -^sup)- Because of this choice, 
it is possible to force resolution down to sublevel I = 4 in- 
side clouds (while the adopted E value ensures subdivision only 
down to / - 3). In total, there are about 1650 clumpfl 

The case of a clumpy disk heated by diffuse radiation is 
shown in the right panel of Fig. |5] (dashed line). The intrin- 
sic radiation emitted by the disk and the bulge is that for the 
Sb template of Fig. H] and the bolometric luminosity Lboi - 
7.3 X lO'^L©. When compared to the model in which all dust 
is in diffuse compo nents (dotted l ine; a nalogous to the two 
dust disk model of iPopescu et alj |2000|) . the clumpy model 
show s the usual reduction in the attenuation of s tarlight (see, 
e.g. IWitt & Gordonlll996l: lvfo)si & Dweklll999h . This is re- 
flected by a reduced dust emission in the infrared, for A ^ 
200yum. At longer wavelengths, instead, the difference between 
a clumpy and a smooth distribution is not appreciable. Both dif- 
fuse dust and dust in clouds externally heated by the ISRF (pas- 
sive or quiescent clumps; iPopescu et al.il2000l) contribute to the 
FlR/submm emission. For the cloud optical depth considered 
here, dust in quiescent clumps is not substantially colder than 
the diffuse dust, and the spectral shape of the two components 
is similar. The FlR/submm emission is thus proportional to the 
total dust mass in a galaxy. The result holds even if E (and thus 
the cloud optical depth) is doubled: the diffuse NIR ISRF is still 
able to penetrate a cloud, and the SED at longer wavelengths 
shows no app reciable change. Similar conclusions are drawn by 
iMisselt et al.l (2001), though they show a more pronounced de- 
pendence of the thermal peak on clumping. The difference is 
probably due to the different heating geometry, with a clumpy 
shell heated by an internal spherical distribution of stars. 

Until now, all the stellar components have shared the same 
spectrum. In the next model, the low-mass spectrum of Fig.|4]is 
used for the disk and the bulge. These components have a total 
luminosity 6 x 1O'"L0, with 70% of the radiation emitted by 
the disk. This is the same luminosity the old component has in 
the clumpy model shown before, which is close to the observed 
stellar SED in the optical and NIR. 

^ Executing the dust emission part of TRADING on the clumpy 
model presented in this section requires about 1.6Gb of RAM on a 64- 
bit CPU to store a grid made of 3.2x10^ cells (90% of which are leaf 
cells). The memory load is mostly due to the pointer-based structure of 
the binary tree. As an increase in resolution results both in an increase 
of memory requirements and calculation time, TRADING efficiency 
will certainly benefit from parallelization with a distributed memory 
scheme. 

^ Simulations with the clumpy dust grid shown here take about 65 
hours on a 2.2 GHz AMD Athlon 64 X2 Dual Core Processor 4200+; 
of these, 58 hours are needed by the dust emission calculation. Maps 
of 512x512 pixels, covering all the dust distribution, are produced for 
three inclinations; A'pho, = 4 x 10^ is sufficient for an adequate S/N of 
the maps. Self-absorption is included and iteration is stopped if the SED 
converge to within 5%. For all the disk models presented here, only two 
iterations were necessary, with the SED actually converging to within 
a few tenths of percent. Self-absorption attenuates dust emission in the 
SED for A S lOyum, especially for the models seen edge-on. However, 
it is found to have little influence on the global energy balance. 



Two UV emitting components are then included, sharing the 
high-mass spectrum of Fig. |4] The first component is a diffuse 
disk, with the same radial and vertical scalelengths as the clump 
distribution. After being extinguished by the adopted dust ge- 
ometry, the spectrum becomes flat, at least for the edge-on case. 
If such component, with luminosity 0.8 x lO'^'L©, is added to 
the old stars, the stellar SED is well matched also at UV wave- 
lengths (Fig.|5] right panel, dot-dashed line). The dust emission 
SED is close to that of the clumpy model with no separate emis- 
sion from high-mass stars. 

A second UV component is represented by point sources 
emitting radiation from the center of each cloud. This compo- 
nent of embedded sources is heavily extinguished; stellar radia- 
tion at small wavelengths is almost completely absorbed and the 
output spectrum peaks at 1 - 2jjm. Because of the high extinc- 
tion, a considerable luminosity can be assigned to this compo- 
nent without modifying the stellar SED of the simulation. The 
luminosity is instead constrained by dust emission: the peak of 
thermal emission can be reproduced with 1.0 x 10^'^ Lq. The fi- 
nal SED is shown with a solid line in Fig.|5](right panel). There 
is a good accordance with observations along all the spectrum, 
with the exception of a lower MIR flux in the simulation. I will 
comment on this later when discussing the dependence of the 
model on the chosen resolution. No significant changes in the 
SED can be obtained if a fraction of the cloud is taken to be 
only passive (e.g. without sources inside); if the UV source is 
positioned randomly in the cloud volume; if the UV luminos- 
ity emitted in a cloud is shared by a few sources. If instead the 
radiation is emitted uniformly inside the cloud, the SED peaks 
at colder wavelengths (an effect similar to a lack of resolution; 
see later). If S is doubled, each cloud becomes smaller and more 
opaque: the SED does not change significantly, with deviations 
within 5% at the peak of the SED and 10% at about lO/jm. 

In total, the model has Lboi = 7.8 x lO'^L©; The UV- 
dominated spectrum of high-mass stars emits about 23% of Lboi 
(sligthly more than in the Sb template), and 55% of this radia- 
tion is embedded inside clumps. By using standard calibrations 
(iKennicuttll 1 99 8*). the intrinsic, dust-free, luminosity at 0.15 fim 
can be converted to a SFRw 3M(T) yr^', ajvalue similar to the UV 
component used in the model of Popescu et al.l ([2000). Of all the 
bolometric radiation, 34% is absorbed globally; 20% is absorbed 
at /l ^ 0.35//m. Despite radiation from high-mass stars suffers a 
larger internal extinction, still radiation from low-mass, older, 
stars dominat e the dust heating. This is instead in contrast with 
Pop escu et al.L where UV radiation dominates the energy bal- 
ance, contributing to 70% of the FIR radiation. The reason for 
this discrepancy is unknown. Since the derived SFRs are sim- 
ilar, it should be due to model differences in the optical/NIR. 
In the current model, extinction of optical radiation could have 
been overestimated at around A = 0.5;um, since the stellar scale- 
length is taken to be about 3 0% smaller than wh at deduced from 
B and V-band observations ( Xilouris et al.lll99^ ; this makes the 
stellar disk more concentrated with respect to the diffuse dust. 
However, most of observed radiation comes at longer wave- 
lengths, where the adopted value for /is is closer t o the fi ts. 
Furthermore, the second dust disk of IPopescu et al.l (l2000h is 
likely to absorb more radiation than the clumpy component in 
this work. 

The cuiTent model also differs i n the fraction of UV radia- 
tion which is embedded in clouds. In lPopescu et al.l (l2000h . dust 
emission caused by absorption of UV photons from localised 
sources is modelled using a SED template based on observations 
of a Galactic H ii region. This template consists in a modified 
blackbody with emissivity oc A'^ and T ^ 35 K; it does not in- 
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Fig. 6. Grid resolution convergence check for the emission from 
clouds in the clumpy model. Only dust emission is shown (see 
text for details). 



elude, in the words of the axiihox^, potentially cold emission com- 
ponents that might be expected f mm "parent" molecular clouds 
in juxtaposition to their "offspring " H ii regions. From the point 
of view of the current work, they supply a potentially cold emis- 
sion component by adding the second dust disk. However, in 
iPopescu et al.l (l2000h the calculation for localised emission is not 
done self-consistently with that of the second, diffuse, disk. This 
might have altered the balance of the UV components, and be 
the cause of the discrepancy with the results shown here. 

iDasvra et all ( |2005|) use the IXilouris et all ( Il999h disk and 
increase the dust emissivity, without adding a second disk. As 
a result, the dust temperature reduces and the thermal emis- 
sion SED shifts toward longer wavelengths, thus matching the 
FIR/submm observations. The missing flux at lower wavelengths 
is provided by absorpti on of U V starlight with luminosity (SFR) 
larger than in Popescu ~et al.l (|2000). Indeed, observations and 
models suggest that the dust emissivity could be lar ger, at least 
in de nser environment (see, e.g., the discussion in lAlton et al] 
12004 . However, the larger dust emissivity they derive may be 
simply due to the lack of cold dust emission in the H ii region 
template they use (which is the same as in Popescu et al. 2000). 

For the clumpy grid shown here and with the current ver- 
sion of TRADING, it is not possible to check for the conver- 
gence of the results by simply increasing the resolution within 
clouds. However, a test can be made by comparing: a simulation 
in which only the embedded UV radiation is emitted from the 
center of the clouds; a simulation of a single cloud with the high 
resolution grid used for the benchmark cases (Sect. 12.5b . scaled 
to the total number of clouds. This is shown if Fig. |6] (dust emis- 
sion only). If the cloud resolution is limited to sublevel I = 3 
(dotted line), the peak emission is at larger wavelengths than 
in the high resolution calculation (solid line). The situation im- 
proves by forcing resolution to sublevel / = 4 (dashed line; the 
value used in this paper), since larger thermal equilibrium tem- 
peratures can be sampled by a higher spatial resolution closer 
to the embedded source. Because of the high extinction in the 
UV, the central source does not heat significantly the external 
parts of the cloud. In full simulation these are predominantly 
heated by the diffuse ISRF, for which a coarser cloud and smooth 
medium resolution is sufficient. Thus, the differences around the 
peak emission are smaller than in Fig.|6] This was confirmed by 



tests on models with less memory requirements (smaller number 
of clouds). For the continuum emission at shorter wavelengths, 
however, differences are higher and could explain the low flux 
predicted at 24/im in the SEDs of Fig. |5] 

Instead, the PAH features show a negligible change with 
resolution, since they depend only linearly, for the temperature 
range considered here, on the intensity of the ISRF (see, e.g.. 
Fig. [3] and Draine & Li 2007). The observed flux at 12yum in 
Fig. |5] could be matched by increasing the amount of embedded 
radiation of about a factor two; however, this would also shift 
the peak of the simulated SED at shorter wavelengths. Also, the 
flux cannot be increased by increasing the amount of mass in 
clouds (and thus their number), since the mass is constrained by 
the FIR/submm observations. Thus, the low predicted flux in the 
PAH features may point towards different ISRF intensities than 
what can be obtained here with a simple cloud description. 

The assumption of spherical, homogeneous, clouds is clearly 
a poor description for the complex st ructure of the molecula r 
medium, which appears to be fractal dBlitz & Wimamslll999l) . 
Unfortunately, the current, non-parallel, version of TRADING 
and the limitations in the available computer memory do not al- 
low a detailed study of the effects of such structure on the models 
output. Limited tests have been carried out by assuming that only 
a fraction of the cells within the spherical volume of each cloud 
is occupied by dust associated with molecular gas, while the re- 
maining cells have the same properties of the underlying smooth 
disk. I have run two cases with filling factors // = 0.5 and 0.2. In 
the first case, clouds maintain a structural integrity, being made 
by many interconnected cell s, while in the second the structure 
is looser (iMisselt et al.ll200ll) . The density of the dust associated 
with molecules increases by a factor 1 /// with respect to the 
standard models, and this akeady may pose resolution problems, 
expecially for the // - 0.2 case if the central source happen to 
be within a higher density cell. As expected, more UV radiation 
is able to escape clouds in these cases (though the // = 0.5 case 
is still consistent with observations while the // = 0.2 is only 
higher by about 30%). As a consequence of that (and probably of 
the limited resolution) the flux between 20 and 80 /im is lower 
than in the standard model, by 20% and 40% at 30/im for the 
// - 0.5 and 0.2 cases, respectively. However, there are no sig- 
nificant differences in the optical SED (because emission comes 
prevalently from the diffuse medium), in the PAHs range and in 
the thermal peak at larger wavelengths. 

4.3. Attenuation vs inclination and images 

The global SED shown before have been obtained by integrat- 
ing over simulated images for the edge-on case. However, im- 
ages can be produced at any inclination. In Fig.|7] true color im- 
ages are shown for the clumpy model with embedded sources, 
seen at / = 86°. The upper panel shows an image in the opti- 
cal. Reddening is evident in the extinction lane and is mostly 
due to the diffuse disk, while clumps, due to their high opti- 
cal depth, simply extinguish all radiation in their background, 
regardless of the wavelength. The lower panel shows an image 
in the FIR/submm around the peak of dust emission. Dust in 
clumps is hotter (bluer) in their centers, due to the embedded 
sources, and colder (redder) close to their surface, where the 
heating is mostly due to the diffuse ISRF. Emission from dif- 
fuse dust shows a gradient, which is strong only in the center of 
the disk, where the dust temperature is larger (mostly as an effect 

of the stronger ISRF in the bulge region). 

The two du st disks model of iPopescu et al.l (12000 ; see also 
iTuffs et al ] |2004) has been successful in reproducing the depen- 
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Fig. 7. RGB true color images of the clumpy model, seen from / - 86°. The upper panel shows an image in the optical (Red = 
K-band, Green= V-band, Blue= B-band). The lower panel shows an image in the FlR/submm, around the peak of dust emission 
(Red = 850 yum, Green= 250 yum, Blue= 70 fim). Each panel has an extent of 32 kpc x 3.5 kpc. 
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Fig. 8. Attenuation as a function of inclination. AM is the differ- 
ence, in magnitude, between the attenuated flux and the flux that 
would be emitted in the absence of dust. 

dence of attenuation on inclination, as obtained from the analysis 
of the luminosity function of disk galaxie s with different axis ra- 
tio in the Millennium Galaxy Catalogue dPriver et al. I l2007h . In 
Fig. [8j 1 show the difference (in magnitude) between the total 
magnitude and the magnitude in the no-dust case (the latter be- 
ing independent of ; for the isotropic emission assumed here). 
In the B-band, attenuation for the clumpy model is reduced, as 
expected, with respect to the two disks modefl However, the de- 
pendence on / is rather similar Thus, the clumpy model could in 
principle be consistent with observations, though a more detailed 
study is needed, with the separate analysis of the attenuation of 
the stellar bulge and disk. The reduction in attenuation (at least 
in the B-band) does not seem large enough to substantially al- 

' I note here that the two disks model shown in this section and in 
Fig. \5\ though analogous in concept, is not the same as the model of 
iPope scu et al. (2000). In this paper, it is used only to show the changes 
when passing from a diffuse distribution to a clumpy distribution for the 
dust component associated to the molecular gas. The ma in difference is 
in the radial scalelength adopted for this component. In lPopescu et alj 
( 12000 ) the scalelength is larger, resulting, for a similar dust mass, in a 
disk with a face-on optical depth roughly half of that for the molecular 
component in this work. 




Fig. 9. K-band images for the two dust disks (top) and clumpy 
(bottom) models, seen edge-on. The extent of the images along 
the major axis is 32 kpc. 



ter the conclusions drawn by [Driver et"aL I (120081) on the energy 
balance of the Cosmic SED. 

In the K-band, the attenuation is rather similar both in 
amount and dependence on /. Indeed, when seen edge-on, both 
models show a dust lane of similar depth (Fig.|9]l. In the two dust 
disks model, most of extinction is due to the more centrally con- 
centrated second dust disk, which has a higher opacity (t^"' k 4). 
In the clumpy model, it is due to clouds, which still have a non 
negligible optical depth (t w 1) and, acting more like foreground 
dust, are more effici ent in extinguishing radiation. Using the pro- 
cedure of Bianchi (2007h . the clumpy image could be fitted with 
a smooth disk with t^°- ^ 0.4 (the fitted radial scalelength is in- 
termediate between that of clumps and of the diffuse disk, while 
the vertical scalelength is retrieved correctly). This is at odds 
with observations, as NGC 891 and other edge-o ns do not show 
a pronounced extiiiction lane and have T y°- $ 0.1 dXilouris et alJ 
ll999l:lDasvra et al.ll2005l:[Bianchill2007 '). A fit to a simulated im- 
age in the V-band, instead, retrieves the parameter of the diffuse 
disk: at optical wavelengths this structure is more evident than 
in the NIR and coherent while the clumps behaves as noisy spots 
that increase the fit residuals (analogous to the inhomogeneities 
in real images). In optical wavelengths, thus, it is possible to 
hide dust in clumps, though the underestimation depends on the 
model assumptions (iBianchi et al.l l2000bl: iMisiriotis & Bianchil 
20021). 

In Fig. [To] I show the major axis surface brightness pro- 
files of dust emission for the two disks and the clumpy mod- 
els seen edge-on. Images have been smoothed with a PSF of 
FWHM=16", a resolution similar to what can be achieved, cur- 
rently or in the next years, at 70/im using MIPS aboard the satel- 
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Fig. 10. Major axis surface brightness profiles for the two disks 
(top) and clumpy (bottom) models, seen edge-on. The pro- 
files have been convolved with a gaussian PSF (FWHM=16"); 
smoothed by folding the left and right part of the image (for the 
clumpy model); arbitrarily normalised at 8kpc. As an aid to the 
eye, the dotted segments show the slopes of exponential fall-offs 
with scalelengths = 4kpc and /is/2. 



lite Spitzer dRieke et al.ll2004. at 250/i m using SPIRE aboard 
the satellite Herschel ( Griffin et alj2007h and at 8 50om using the 
submm c ameras SCUBA2 ([Holland et al.ll2006') and LABOCA 
dKrevsa et al. 2003). To these, I have added a profile at about 
8yum, centered on a large PAH feature. For a diffuse dust disk co- 
spatial with the heating sources, as is the case in the two dust disk 
model for the second dust disk and the stellar disk, the emission 
gradient becomes less steep as the wavelength increases. In the 
submm, the Rayleigh- Jeans spectrum does not depend strongly 
on the dust temperature, and the gradient approaches the intrin- 
sic slope of the dust distribution (close to the radial scalelength 
/iFl). Instead, at wavelengths shorter than the peak, the spec- 
trum depends more on the strength of the ISRF, as shown by the 
IQfim profile (the further steepening in the center being due to 
the bulge). This is not the case for the PAH emission, however, 
because of the linear dependence of the spectrum on the ISRF 
intensity. 

When the dust mass of the second disk is distributed in 
clumps (retaining the same radial scalelength), the differences 
between the profiles at aU wavelengths are reduced. This is 
caused mainly by the contribution of internal clump emission 
at A < 200yum, which is independent of the clump position. 
Also, the Sjum profile becomes closer to the submm, explain- 
ing the observed cor relation between PAHs and cold dust (see, 
e.g. lHaas et alj20()2h . Though reduced, the bulge contribution to 
the lOfim emission is still detectable (see also Fig.|7J. 

Finally, I show in Fig.[TT]a comparison between the clumpy 
model and the SCUBA major axis profile as presented in 



The radial profile at 850yum is steeper than the exponential of scale- 
length As, both because of the contribution of the diffuse disk, and 
because the projected profile of an exponential disk seen edge-on is 
oc rKi(r/h), with Ki{x) the modified Bessel function of the second kind 
and first order and r the distance from the center ( Kylafis & Bahcall 
1 19871') . For an infinite disk, the profile tends to the exponential. In the 
case shown here, the second dust disk, the distribution of clumps, and 
the stellar disk are all truncated at 4 As. This causes the profile steepen- 
ing seen at 16 kpc. The truncation has little effect on the results. 
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Fig. 11. Major axis surface brightness profile at 850 yum for 
the clumpy model , compared to the SCUBA observations of 
lAlton et al.l ( |2000l) . The model has been convolved to the 
SCUBA beam (FWHM=16") and the profile smoothed as in 
Alton et al. (2000). 



lAlton et alj (|2000|) . As for model of'Popescu et al.' ^2000^ (see 
their Fig. 6) there is a broad agreement with observations, both 
in flux level and gradient. Only, the clumpy model appears to 
be smoother than observations, possibly hinting to a more com- 
plex structure in the cloud distribution than what adopted for this 
paper. 

5. Summary 

I have presented in this paper the Monte Carlo radiative trans- 
fer code TRADING. It is designed to study continuum extinc- 
tion and emission in a dusty medium. The main features of 
TRADING are: 

i. a description of the dust distribution based on a binary-tree 
adaptive grid; 

ii. possibility of emitting stellar radiation from several diffuse 
distributions and collections of point sources/clouds, each 
with an independent spectrum; 

iii. dust emission from a distribution of grains of different sizes 
and materials, under thermal equilibrium or stochastic heat- 
ing conditions; 

iv. self absorption through an iterative process. 

The code shows an excellent agreement with benchmark 
cases and can be easily adapted to several environments of as- 
trophysical interest. For the application shown in this paper, 
TRADING has been tailored to the case of a galactic disk. The 
disk includes: a diffuse dust component; a clumpy component 
made of a distribution of spherical clouds, with properties mod- 
ell ed on those of GMC in loc al galaxies. Bu i lding on the studies 
of IXilouris et alj fl999) and iPopescu et al.l (l2000l) . simulations 
have been produced for the well known case of the edge-on spi- 
ral NGC 891. The main results can be summarised as follows. 

1 . The global SED of NGC 89 1 from the UV to the mm can be 
reproduced by a model including: a diffuse dust disk, more 
widespread than the stellar and of moderate optical depth 
(Ty° =l); a clumpy dust distribution, of mass similar to that 
of the diffuse disk; a population of old stars, smoothly dis- 
tributed in an exponential disk and a bulge, emitting about 
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75% of the bolometric radiation; a population of young stars 
emitting the remaining fraction mainly in the UV, roughly 
half of which embedded in clouds. Despite the increased at- 
tenuation for UV radiation, the old stars are the main con- 
tributor to dust heating (about 60% of all ra diation emitted 
by dus t). This is in contrast with the result of iPopescu et al.l 
(l2000h . which predict a dominant UV contribution to heat- 
ing. 

In the MIR, the current model produces less radiation than 
the observations. This can be only in part explained with a 
lack of spatial resolution for clouds. The difference can be 
possibly due to a more complex structure of the ISM and 
ISRF than modelled here. This could also be the reason for 
the larger extinction in the K-band image and the smoother 
appearence of the 850//m than observed. 

2. Dust in clouds is found to be significantly heated by both 
the embedded UV radiation and the NIR radiation from the 
diffuse stellar distribution penetrating the outer shells of a 
cloud. This second source of heating results in cold dust 
emission which is not present in the template for localised 
emission used by f opescu et al. (2000) and iDasvra et al.l 
(l2005b . In iPopescuet al.l (12000). cold dust emission is pro- 
vided by a second diffuse dust disk; however, contrary to 
this work, the radi ative transfer c alcu lation for the extra 
dust component in iPopescu et al.l (12000) is not done self- 
consistently with that for locahsed sources. This is the cause 
of the larger fraction of embedded UV radiation in the cur- 
rent work. Also, the difference between a clumpy and smooth 
disk results in a smaller attenuation of diffuse radiation at 
optica l wavelengths in the current model. In Dasvra et alj 
(I2005h . the lack of cold dust emission associated with lo- 
calised emission appears to be the main reason of their claim 
for a larger dust emissivity in the FIR. Other models which 
compute the contribution of clouds indep endently from tha t 
of the diffuse ISM, such as GRASIL dSilva et al.l Il998l) . 
might need to include external heating of clouds by a dif- 
fuse ISRF. Because of the diffuse heating, cold emission in 
the submm comes from all the dust components. 

3. In a clumpy disk, and for the resolution of current and future 
instruments, color gradients are rather similar along all the 
dust emission spectrum. They follow the gradients of stars 
and of the main dust component. A steepening may be ob- 
served in the Wien side of the spectrum for thermal equilib- 
rium heated grains (20jum ^ A ^ 100/zm), mostly because of 
the presence of a bulge. 

The conclusions drawn on the energy budget for the single 
object NGC 891 might not be of general application for other 
galaxies, though NGC 891 is considered to be a prototypical 
spiral. Unfortunately, not many galaxies have detailed obser- 
vations of the FIR/submm beyond the thermal peak. However, 
the advent of the SP IRE instrument aboard the Herschel satellite 
(iGriffin et al .1120071) . in conjunction with the new ground-based 
submm facilities, will soon permit a more detailed study of this 
spectral range, necessary for constraining the bulk of dust mass. 
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